#ifndef _DWMRICORE_DTI_H
#define _DWMRICORE_DTI_H


#include "Point3F.h"
#include "defines.h"

#include <stdlib.h>

// TEEM
#include <teem/nrrd.h>
#include <teem/biff.h>


/********************************************/
/* data in m_pDTI: flag xx xy xz yy yz zz	*/
/********************************************/
class CDTI
{
public:
	float *m_pDTI;
	uint m_w;
	uint m_h;
	uint m_d;
	uint m_n;			// should always be 7
	uint m_order;		// should always be 2

public:
	CDTI(void);
	~CDTI(void);

public:
	// open the DTI in Nrrd format
	int OpenDTIFile(const char *pathname, FILE *fp = stdout);
	// save the DTI as Nrrd format
	int SaveDTIFile(const char *pathname, FILE *fp = stdout);
		
	int CreateDTI(const uint w, const uint h, const uint d, const float *dti = NULL);

	// return the pointer of the DTI data
	float *GetDTIData() { return m_pDTI; }		

	// use the following interface to get the DTI value 
	// of a particular position in the volume
	void GetDTI(const uint x, const uint y, const uint z, float *dti);
	void GetDTI(const float x, const float y, const float z, float *dti);
	void GetDTI(const CPoint3F &p, float *dti);
	void GetDTI(const CPoint3F *p, float *dti);

	void GetDTIFast(const uint x, const uint y, const uint z, float *dti);
	void GetDTIFast(const uint index, float *dti);

	void SetDTI(const uint x, const uint y, const uint z, float *dti);

	// use the following interface to get the size of the DTI volume
	uint GetWidth() const { return m_w; }
	uint GetHeight() const { return m_h; }
	uint GetDepth() const { return m_d; }
	uint GetN() const { return m_n; }
	uint GetOrder() const { return m_order; }

	bool HasData() { 
		if (m_pDTI) {
			return true; 
		} else {
			return false;
		}
	}

	void DTI2Direction(const float *dti, CVector3F *dir) const;
	float DTI2FA(const float *dti) const;
	float DTI2CL(const float *dti) const;
	float DTI2CP(const float *dti) const;
	float DTI2MD(const float *dti) const;
	void DTI2Eigen(const float *dti, float *eval, float *evec) const;		/* compute eigen value and eigen vectors */
};

#endif
